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Abstract 

A simple algorithm is described to sample permutations of identical particles in Path Integral 
Monte Carlo (PIMC) simulations of continuum many-body systems. The sampling strategy illus- 
trated here is fairly general, and can be easily incorporated in any PIMC implementation based 
on the staging algorithm. Although it is similar in spirit to an existing prescription, it differs from 
it in some key aspects. It allows one to sample permutations efficiently, even if long paths (e.g., 
hundreds, or thousands of slices) are needed. We illustrate its effectiveness by presenting results 
of a PIMC calculation of thermodynamic properties of superfluid 4 He, in which a very simple 
approximation for the high-temperature density matrix was utilized. 

PACS numbers: 02.70.Ss, 05.30.-d., 67.40.-w 
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I. INTRODUCTION 



The Path Integral Monte Carlo method is arguably the most powerful numerical technique 
to calculate thermodynamic properties of continuum (i.e., non-lattice) quantum many-body 
systems at finite temperature.— For Bose systems, it is the only known, generally applicable 
theoretical method essentially free from approximations. Numerical estimates yielded by 
PIMC are affected by a statistical error, as well as by systematic errors, due to the finite 
size of the simulated system and to imaginary time discretization. In most cases, however, 
the computational resources typically available nowadays allow one to render the size of all 
of these errors insignificant in practice. 

The most notable application of PIMC to date, is the study of the superfluid (SF) tran- 
sition in liquid 4 He by Ceperley and Pollock, 2 whose results have become the standard 
reference for all theoretical calculations on SF helium; but numerous applications to other 
Bose systems have been reported in the literature, over the past two decades. 
No general formulation exists as yet of PIMC (nor of any other Quantum Monte Carlo 
method), capable of overcoming the sign problem, that has so far made it impossible to 
obtain equally high quality results for Fermi systems. Even for fermions, however, PIMC 
proves a valid option, allowing one to obtain approximate estimates, of accuracy at least 
comparable to that afforded by other methods.— 

Physical effects of interest in quantum many-body systems are almost invariably asso- 
ciated with quantum statistics; for example, superfluidity in 4 He is intimately connected 
to long exchange cycles of helium atoms. Because a direct summation of all N\ permuta- 
tions of N indistinguishable particles is unfeasible, except for very small values of N, within 
PIMC quantum statistics is included by performing a statistical sampling of permutations. 
Thus, an all-important ingredient of any practical implementation of PIMC is an efficient 
procedure to carry out such a sampling. 

Since the pioneering study of Ref. Q, there has been relatively little experimentation 
with implementations of PIMC differing, in some of the more important aspects, from the 
one described in Ref. ^, henceforth referred to as CP. The CP implementation has come 
to be regarded as "canonical", especially when studying quantum many-body systems in 
the highly degenerate regime (i.e. at low temperature). It is based on an accurate ("pair- 
product") high-temperature density matrix, allowing one to observe convergence of the 



2 



physical estimates with a relatively low number of imaginary time "slices" (of the order 
of 40 for superffuid 4 He at a temperature T—l K). Two slightly different procedures have 
been proposed and utilized, in the context of CP, to perform the sampling of permutations, 

n 

both of which are thoroughly described in Ref. UJ. To our knowledge, no systematic, quan- 
titative assessment of the relative merits and advantages of these two sampling strategies 
has yet been offered; it is also unclear to what extent their effectiveness and applicability 
are problem-dependent, and/or hinge on the use of the above-mentioned high-temperature 
density matrix. 

In this work, we illustrate a new method to sample permutations of indistinguishable 
particles in PIMC simulations. It bears some similarities with one of the two strategies 
described in Ref. Ill, but differs from it in some important technical aspects. We also deem 
it easier to implement, and may be potentially more efficient, even though, naturally, this 
speculation will need to be supported by systematic comparisons with the other existing 
options. 

As an illustrative application of our sampling method, we have carried out a PIMC 
simulation of liquid 4 He in the SF regime, i.e., we have repeated the original calculation of 
Ref. 2). SF helium is the accepted test bench for quantum Monte Carlo calculations, since 
it is the most extensively studied quantum fluid, for which effects of quantum statistics 
manifest themselves at the macroscopic level. In order to make the test more significant, 
and help expose any deficiency or merit of the permutation sampling procedure, we have 
not utilized the pair-product high-temperature density matrix; rather we opted for a much 
simpler form, which requires a substantially larger number of imaginary time slices, in order 
to observe convergence of the estimates. Besides providing results for energetic properties, 
known to be affected quantitatively by Bose statistics, and which we compare to those of 
Ref. 12, we have also attempted to furnish here some quantitative information, which should 
help assess the efficiency of our method in sampling the space of all possible entangled many- 
particle paths (i.e., including permutations). It is our hope that this will provide a baseline 
for future, more extensive comparisons of different approaches. Somewhat interestingly, our 
results show that the PIMC simulation of SF 4 He is feasible, albeit at a higher computational 
cost, even with a relatively simple PIMC implementation. Doubtless, this is also in part 
due to advances in computing hardware, which enable what may have been prohibitive two 
decades ago, when the first such simulation was carried out. 
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The remainder of this manuscript is organized as follows: in the next section, we provide 
a detailed description of our computational methodology. In the following sections, we 
illustrate our results, and outline our conclusions and outlook. 



II. METHODOLOGY 

The PIMC methodology is fairly mature, and extensively described in Ref. Q, to which 
we refer readers interested in a thorough, comprehensive illustration. Our specific imple- 
mentation is largely based on the ideas and methods presented therein. Nevertheless, at the 
risk of some redundancy, we provide a somewhat detailed description of our implementation 
here. This will hopefully facilitate the task of others who may wish to repeat our study 
and/or experiment with our algorithm to sample permutations. 

Consider a quantum many-body system of N identical, point-like particles of mass m, 
described by the following Hamiltonian: 

JV 

^ = -AEV l 2 + En|r,-r,|) (1) 

8=1 i<j 

where A = h 2 /2m. Implicit in the above model is the assumption that the interactions 
among particles can be accurately represented by a pairwise, central potential (the V term 
in (P)). Although this is obviously an approximation, it is commonly made in theoretical 
studies of most quantum fluids. In any case, it is not a requirement for the applicability of 
PIMC, nor of our specific implementation. In the following, it is assumed that particles in 
the system obey Bose statistics.- The system is assumed to be enclosed in a vessel, shaped 
as a parallelepiped, with periodic boundary conditions in all directions, and to be held in 
thermal equilibrium at a temperature T. 

The thermodynamic average of a physical quantity formally represented by an operator 
O (for simplicity assumed diagonal in the coordinate representation) is expressed as follows: 

(0) = i J dRO(R) p(R,R,(3) (2) 

where (3=1 /T (we work with units where the Boltzmann constant ks—l), and R = 
{r%, r 2 , ...rjv}, is a configuration of the system, specified by the positions of all the N 
particles. In Eq. (j2J), p(R,R',(3) = (_R|e _/3H |-R') is the many-body density matrix, and 
Z = J dR p(R, R, (3) is the partition function. 



A Monte Carlo evaluation of (J2J) consists of generating a large set of random many- 
particle configurations {R p }, p = 1,...,M, statistically sampled from a probability density 
proportional to p(R,R,f3); the thermal average (J2J) can thus be estimated as a statistical 
average over the set of values {0(R P )}. 

An explicit expression for the density matrix p(R, R, (3) is unavailable for any non-trivial 
many-body system; however, one can still generate the set {-Rp}, by sampling discrete many- 
particle paths X p through configuration space, i.e., 

X p = {Ro p , Rip, R2 P , ••• Rbp] (3) 

Paths are formally defined in the imaginary time interval < r < (3, i.e., Rj = R(j5r), 
with L5t = f3, and are randomly drawn from a probability distribution p(X) given by 

L-l 

p{X) = p(Ro, R u R 2 , ...R L ) = J] Po(Rj, R j+ i, St) (4) 

3=0 

where p Q is an (analytically known) approximation to the true many-body density matrix, 
constructed to be asymptotically exact in the "high temperature" St — > limit. It can be 
shown that in that limit {L — > oo), each configuration R p visited by paths is statistically 
sampled from a distribution proportional to p(R,R,(3). 

The configuration R(f3), i.e., that corresponding to the end of the path in imaginary 
time, must coincide with R(0), except for a permutation P of the particle labels (1 through 
N) . The possibility of permutations of particles must be allowed, in order to incorporate in 
the calculation the effects of particle indistinguishability and Bose statistics. Consequently, 
although many-particle paths are periodic in imaginary time, i.e., the configuration R(J6t+ 
q(3) (q being an arbitrary integer) is identical with R{J8t) (in that particles occupy identical 
positions), individual particle labels can be different. Stated differently, single particle paths 
t i(JSt)---Tn(JSt) can become "entangled", as a result of permutations. 

Permutations normally become important at sufficiently low temperature; at high tem- 
perature, only the identity permutation contributes significantly to thermal averages. At 
low temperature, however, permutations underlie phenomena such as superfluidity and Bose 
Condensation in liquid helium and, presumably, in all other superfluids. 

In an actual calculation implementing the above computational scheme, one must neces- 
sarily work with a finite value of L; in principle, one ought to carry out the L — > oo limit 
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by extrapolating numerical results obtained with different values of L. In practice, this pro- 
cedure proves quite cumbersome, especially when one is interested in many thermodynamic 
points. Thus, one typically performs all calculations (at a given temperature) with a single 
value of L, chosen sufficiently large so that estimates may be expected to coincide with the 
extrapolated values, within some small tolerance. For reasons of efficiency, it is desirable 
that such "optimal" value of L not be too large (a few hundred slices at the most); thus, it 
is advantageous to work with as accurate as possible a "high-temperature density matrix" 
p , which will allow one to achieve convergence of the numerical estimates without resorting 
to impractically large values of L. 

The importance of this issue was demonstrated by Ceperley and Pollock,— who proposed the 
following form for p a : 



p a (R, R', 6t) = A F (R, R\ 6t) I J] exp -u(r y , r' 5r 



i<j 



(5) 



where 



Tj — r-j, A F (R,R',5r) = ]lj=i PF( r i, r i, ^ T ), is the exact density matrix of a 



system of N distinguishable, non-interacting particles, with 



-3/2 



exp 



IV — r 



t\2 



4X5r 



(6) 



and where u is obtained by imposing that p a be the exact density matrix for a system of 
two interacting particles. For PIMC calculations of highly quantal, hard-sphere-like systems 
such as condensed Helium, the form (j3J) for p Q affords a tremendous increase in efficiency, 
with respect to other, simpler forms for p Q (such as the so-called primitive approximation; 
for details, see Ref. Q). 

In this work, we have not made use of the high-temperature density matrix (JSJ), choosing 
instead the following form: 



p (Rj,R j+ i,5r) = Ap{Rj, Rj+i, St) exp 



where 



U(R S ) 



2V(R J ) 



+ V{Rj 



(7) 



(8) 



V(R) = J2i<j V(\ r i — r j\) being the total potential energy of the system in the configuration 
Rj, and 



(9) 



1=1 



if j is odd, and zero if j is even. Here, ViV(R) is the gradient of the total potential energy for 
the configuration R, with respect to the coordinates of the zth particle. This is a particular 
case of a more general expression, which can be shown^^ to be accurate up to terms of order 
t 4 in the expansion of the exact density matrix p(R, R', r) in powers of r. 

Using the form (jSJ) instead of the (superior) pair-product approximation (j5J), results in 
a substantially larger value of L required to achieve convergence. For the specific physical 
system that we have chosen to test our algorithm, namely superfluid 4 He, the number L 
of imaginary time slices needed is as much as 16 times greater than if (0) had been used. 
The reason for our choice is that our interest in primarily methodological. Specifically, we 
wish to separate the relative contributions to the effectiveness of a PIMC implementation, 
of the permutation sampling procedure and of the high-temperature density matrix utilized. 
A more stringent test is provided of our sampling scheme, if it can be shown to work 
satisfactorily with a fairly simple approximation for p Q . 

A. Path Sampling 

The generation of the set of many-particle paths {X p }, with p = 1, 2, ...M, can be conve- 
niently carried out using the Metropolis algorithm. According to the standard procedure — 
one performs a random walk through the space of A-particle paths A, defined above, start- 
ing from an initial point X . The A p 's are then the points sequentially visited by the random 
walk. 

Let Xi be the Ith element of the set {A p }; in order to generate A^ +1 , one samples a 
modification of the path Xi, involving new positions of one or more particles at several 
points (i.e., configurations) along A/. Let Xf be the path arising from such modification 
of Xi, and let T{X\ — > Xf) be the probability with which Xf is sampled from X\. The 
proposed new path is accepted, thereby becoming the next point of the random walk (as 
well as the next element X\ + \ of the set {X p }), with probability 

W{Xl Xi) = ~pVQ t{ Xi - x t ) (10) 

This is simply done by drawing a random number \ between zero and one; if W[X\ — > 
Xi) > x, then X t+1 = Xf, otherwise X i+1 = X x . 

Of fundamental importance to the efficiency, unbiasedness and correctness of the algo- 
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rithm, are the elementary moves whereby one generates the "trial" path X* starting from 
X\. In our PIMC implementation, two different types of moves are performed. A detailed 
descriptions of these moves is offered in the next two subsections. 



1. 



'Wiggle" type moves 



These moves modify the current path X\ by just altering the path of one particle, ran- 
domly chosen. Random displacements are applied to a number s — 1 of consecutive positions 
of that particle along its path. This can be thought of as "chopping off" a portion of path, 
and replacing it with a different segment. The maximum number of positions modified by 
the update is L — 1, as the two ends of the portion are left unchanged. Because paths are 
periodic, it is possible to update a portion of path of a single particle that will include the 
zeroth or Lth positions. 8 

In order to illustrate this type of move in detail, let us assume that a particle has been 
selected, and let the portion of path to be updated include the positions r fc+1 ...r fc+s _ 1 , where 



< k < L — lisan integer number, and s = 2 m is chosen so that 2 < s < L. Let r' h 



l k+s-l 



be the tentative new positions of the particle, selected according to some (yet unspecified) 
probabilistic criterion, expressed by a sampling function T. For notation purposes, we also 
define r' k = r k and r' k+s = r k+s . Based on (JHJ) and (fTUj) . the acceptance probability of the 
move will be 



W 



{lljJo pF(*' k+j , r'k+j+i, St) 


• exp 




T(X' - 


■+X) 


|nj=o PF(r k+j ,r k+j+1 ,5T) 


> exp 




T(X - 


+ X') 



having defined R' as the configuration that differs from R only by the displacement of 
the chosen particle from r to r', X = Rq, Ri, ...Rl is the current path, whereas X' = 
R , Ri, ...R k , R' k+1 , ■■■R' k+S _i, Rk+s, ■■■Rl is the proposed new path. There is considerable 
freedom in choosing the sampling probability T,— but it is clearly advantageous to do so in 
a way that will simplify the expression (fTTj) . An obvious choice is 



T(X -> X' 

which reduces the acceptance probability (fTTj) to 

W = exp 



s-1 
3=0 



St) 



-5r S "£(u(R' k+J ) 
i=i v 



U(R k , 



(12) 



(13) 
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The probability T so defined can be conveniently sampled by means of the "staging" 
algorithm. 9 The idea is as follows: T is given by a product of s = 2 m Gaussian terms, 
namely 

\2- 



T(X -> X') oc exp 



k+ls 



AX5r 



x exp 
. x exp 



( r fc+l r 'k+2) 2 



4X5 



T 



fe+s-1 



\2 n 



AX5r 



(14) 



Using some simple algebra^ it is possible to re-organize this product in the following, 
"hyerarchical" form: 



T oc exp 



-k+s) 



\2i 



x < exp 



4sX5t 
hi ^2 



x exp 



1 k+s/2 



r 'k,k+s) 



sX5i 



-fc+s/4 L k,k+s/2) 



s\8t/2 



x < exp 



fc+s/8 1 fc,fc+s/4> 1 



exp 



-fc+5s/ 



sA5r/4 

: fc+s/2,fc+3s/4) 



sA5r/4 



exp 



exp 



exp 



l fc+3s/4 1 fe+s/2,fc+s/ 



sA5r/2 



L fc+3s/ 



— r 



fc+s/4,fc+s/2^ 



sA5r/4 

C k+7s/8 ~ r k+3s/4,k+s) 



\2, 



x < exp 



-fc+s/16 



sA5r/4 

^2 



... etc ... 



(15) 



where we have defined r' 



(r' a + r'p)/2. All distances are assumed to be computed 
with periodic boundary conditions, using the minimum image convention. Expression (|15|) 
immediately suggests a sequential, multi-level procedure to generate trial random positions 
of the particle being displaced. Since = and r' k+s = rk+ s , the first factor in (|T5|) 
does not enter the sampling in this type of move. Thus, one starts by generating a new 
"midpoint" position r' k+s ^ 2 , by sampling a three-dimensional Gaussian distribution function 
of semi-width ctq = 



'sXSt/2, centered at (r' k + r' k+s )/2. It is customary to refer to the 
generation of the new midpoint as the zeroth level (1 = 0). 

One then proceeds to the first level (/ = 1), where the two random positions r' k+s ^ 4 



and r' k+3s ^ 4 are sampled from Gaussian distribution functions of semi-width o\ = ysX5r/A, 
centered at positions (r' k + r' fc+s / 2 )/2 and (r' k+s / 2 + r' k+s )/2. At the next level (/ = 2), one 
generates four new positions and so on. The hh level involves the generation of 2 l new 
positions, sampled from Gaussian distribution functions of semi-width 07 = yj2~ l ~ 1 s\8r. 
This "bisection" procedure ends when new positions of the particle have been generated at 
k + 1, k + 2, k + s — 1. The last level is obviously I = m — 1. 



The proposed new path X' may be either accepted or rejected following an acceptance 
test based on (|13p. It proves much more efficient, however, to break down this global, 
final acceptance test, into m — 1 intermediate acceptance tests, each following every level 
of update. Specifically, after completing the Zth level one proceeds to the next level with 
probability 



aborting the the process (i.e., rejecting the proposed new path in its entirety) on the first 
negative outcome of an acceptance test. It is simple to see that the overall acceptance prob- 
ability for the new path remains the same, given by on breaking down the acceptance 
test by levels as explained above. The improvement in efficiency comes from the fact that the 
final acceptance is mostly influenced by the largest displacements, e.g., that of the midpoint. 
Thus, one can reject early (i.e., after the first level), and with relatively little computational 
effort, moves that most likely will eventually be rejected anyway. 

The value of s (namely the length of the portion of path that is updated) is set to ensure an 
optimally efficient sampling. The minimum possible value is s = 2, which gives the highest 
acceptance rate, but at the same time also produces a modest path update (a single point 
of the path is modified). This becomes inefficient at low temperature, as paths can be fairly 
long (e.g., several hundred slices) and such "single-slice" updating can result in a very slow 
diffusion through configuration space, and consequently in undesirably long equilibration 
and auto-correlation times. It is therefore advantageous to take s as large as possible, while 
still ensuring a reasonably high acceptance rate for multi-level moves (acceptance rate is 
a rapidly decreasing function of s). In our algorithm, we typically adjust s so that the 
acceptance rate remains roughly between 20% and 50%. 

2. "Permute" type moves 

These moves involve a group of 1 < n < N particles. They are similar to the wiggle type 
moves, in that corresponding portions of the paths of the n particles are modified, at s — 1 
consecutive points. An additional aspect, however, is that the modified portion of the path 
of a particle in the group will connect, at k + s, to the path of a different particle, among 
the n selected. This elementary move clearly allows one to sample permutations of the N 
particles in the system, over the imaginary time interval [0,/?]. 




j e level I 



E 




(16) 
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The basic scheme of the move is as follows: first, a permutation of particle labels at k + s 
is sampled; the number n of particles involved in this permutation (henceforth referred to 
as the cycle) is not chosen a priori, but can vary from 2 to N. Once the permutation is 
selected, new single-particle paths are constructed, in a way completely analogous to that 
used in the "wiggle" moves. Finally, the new many-particle path X' so obtained is accepted 
with probability given by ([lip. Again, it proves convenient to choose a sampling probability 
for the permutation that will in turn simplify the acceptance probability. Going back to 
Eq. (|15jh the first term of the product is now used to sample permutations, whereas the 
remaining terms are used to construct paths consistent with the permutations that have 
been sampled. 

The sampling of a permutation is a recursive process in which particles are successively 
added to the cycle. The addition of a single particle includes an acceptance test, and the 
sampling of the particle from a table. One begins by selecting a random particle, say the 
z/th for definiteness. Based on (fT5|) . a table is constructed, Kj$, with entries as follows: 

K il} = pF(r uk , r^k+s, s5t) (1 - 5 UU ) (17) 

where r„& is the position of the uth. particle particle at point k, whereas r^k+s is that of 
the a>th particle at point k + s.~— At this point, a first acceptance test is performed, namely 
the process will continue on to the next stage (i.e., selection of the permutation partner for 
particle v) with probability 

(7(1) = 9l (18) 

Qi + p F (r uk ,r uk +s,s5T) 

where Q\ = ^if^. If the acceptance test fails, then the process is aborted, i.e., the 
permutation move is rejected. Suppose, instead, that a positive outcome is obtained; an 
entry a is then sampled from the table with probability U a = K^/Q\. We see 

from (fTTjl that particle v itself is sampled with probability zero, i.e., the sampling of a 
"non-identical" permutation is forced here. The particle labeled a is selected as the second 
member of the permutation cycle being constructed. That means that, in the trial path X', 
the path of particle v will go through r u k at the kth point and through r a k+ s at point k + s. 
At this point, one has to sample a new position of particle a at k + s. Just as for particle 
u, one constructs a table 

K al = PF(r ak , r^k+s, s5t) (1 - 5 auj ) (19) 
11 



and another acceptance test analogous to (|T5j) is carried out, based on the probability 

C(2) = 7\ FT (2°) 

with Q 2 = Z)w^oS (Again, the process is aborted if this acceptance test fails). An entry 
[i is sampled with probability 11^ = K^/Q 2 ; in this case, particle v is sampled with finite 
probability, as the cycle can close on the initial particle, whereas particle a is now excluded 
from the sampling. At this point, if fi — v, then the permutation cycle is closed, and it 
includes two particles, namely v and a. If, on the other hand, \l ^ v, then one must find 
another particle 7, which will become a member of the cycle, such that = r 1 k+s- Again, 
one constructs a table as above, the only difference being that now both a and ji are 
excluded from consideration, as both r afc+s and r Mfc+s are already taken: 

K< $, = PF^^k, r wfc+s , s5t) (1 - 5 auj - 5^) (21) 

A new acceptance/rejection test is performed, based on a probability defined analo- 
gously to CEED-flU: 

C (3) = ^ (22) 

Q 3 + PF(r Alfe , s5t) + piKVc, r ak+s , s5t) 

with Q 3 = J2ui Km} j an d in case of success, one proceeds to sample entry 7 from table K ^ , 
with probability II 7 = K$/Q 3 . 

The basic idea should now be clear: This procedure is iterated until the cycle is finally 
closed, namely until particle v is obtained from the sampling of the table i^ n_1 ). Two 
fundamental aspects of the above scheme to sample permutation cycles, are the exclusion 
from the tables of entries corresponding to particles already in the cycle (the uth is only 
excluded from K^), and the acceptance tests based on C^ n \ preceding each new particle 
selection. The sum at the denominator of C^ 1 ' includes, besides Q n , free-particle density 
matrices associated to all the entries excluded in the sampling table K^ n \ 

Once a complete cycle has been obtained, one must construct trial paths for all particles 
in the cycle, consistent with the selected new positions at slice k+s. This second part is done 
in exactly the same way as for the "wiggle" moves, using the same sampling probabilities. 
Specifically, new midpoint positions are first sampled for all particles in the cycles; then, 
new positions at k + s/2 and k + 3s/4 are sampled, and so on, with acceptance tests as in 
(Unj) after each level of path update. Note that the values of U at points k and k + s remain 
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unchanged, as no particle is displaced at these points; only particle labels are altered at 
point k + s. It is a simple matter to show that the path sampling probability arising from 
the above scheme is indeed consistent with (II (J I) . 

n 

The above scheme to sample permutations is similar to one described in Ref. UJ; the main 
difference, possibly significant, is that, in our procedure, one need not include, in any of the 
acceptance tests, a sum of terms representing all n starting points of the cyclic permutation. 
Moreover, in our method n distinct particles are sampled by construction, and the "identity" 
permutation, namely the one in which particle labels are left unchanged from k to k + s, is 
excluded from the sampling. 

Just as in the "wiggle" moves, s must be chosen appropriately, namely, long enough that 
non-trivial permutation cycles can be sampled with appreciable probability. If s is small, 
particularly if one is working with a small value of 5t, the functions pF are negligible small 
for distances of the order of the average distance between particles, rendering it exceedingly 
unlikely to go beyond the first acceptance test (Eq. IT8|) . On the other hand, taking s too 
long, while allowing for large permutation cycles being sampled, results in very low overall 
acceptance for these cycles, much for the same reasons why acceptance falls for the "wiggle" 
moves as well if s is too large. In the calculations whose results are illustrated in the next 
section, we have generally found that the optimal choice of s is generally the same as for 
the "wiggle" moves. 

Even when s is optimally chosen, typical values of acceptance for permutations are low. 
One tries to keep the efficiency reasonable by attempting a large number permutational 
moves, which can be done fairly rapidly within the relatively simple scheme outlined above. 
We typically attempt several tens of thousands of permutations between two consecutive sets 
of "wiggle" type moves, in which full updates of the paths of all particles are attempted. 

III. RESULTS 

We now describe the results of our PIMC simulations of condensed bulk 4 He at low 
temperature (1 K < T < 4K), obtained with the algorithm described in the previous section. 
The model Hamiltonian for the system of interest is given by |Q. For the purpose of 
comparing our results with existing calculations, 1 ' 2 we used an early version of the Aziz 
potential to describe the interaction between a pair of 4 He atoms (A=6.0596 KA 2 ). 
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We have computed several energetic and structural properties of the system. We have 
observed convergence of the energy estimates with a value of the "imaginary time step" 
<5r=l/640 K" 1 . All of the results presented in this section are obtained with this value of St. 
We estimate any residual, systematic error on the energy arising from our path discretization, 
to be worth no more than 0.1 K per 4 He atom. For structural properties, we have observed 
that estimates obtained with (up to four times) larger values of St are indistinguishable, 
within statistical uncertainties, from those obtained with the above-mentioned value of St. 

We found our optimal value of s, for both "wiggle" and "permute" type moves, to be 
s = 6.— Accordingly, the length of the portion of path that is updated on each move 
is 2 s = 64 imaginary time slices, which corresponds to an imaginary time interval of 0.1 
K _1 , with our choice of time step. Unless otherwise stated, the number of particles in our 
simulated system is iV=64, as in Ref. 0, but results for N=21Q were obtained as well.. 



A. Energetics 

The energy estimators utilized in this work are described, for instance, in Ref. 0. Specif- 
ically, the average kinetic energy per particle K is obtained as 

<*> H Wr ~ lM r * " r -)) + ^(( W ^) 2 ) (23) 
where (...) stands for statistical average, (r& — r^+i) 2 is the square distance between the 
positions of a particle at adjacent points along the path, whereas the gradient of the potential 
energy in the third term is taken with respect to the coordinate of one of the particles, at 
an "even" slice. The potential energy per particle (v) is instead obtained as 

(v) « (24) 

Both relations (|2*Hj) and (J2*3j) are approximate, approaching the exact results only in the 
limit L — > oo, St — > 0. The above kinetic energy estimator is not the most efficient; it 
is known that the so-called "virial" estimator yields more accurate results (i.e., smaller 
statistical errors), given the same amount of computing time.— However, the estimator 
(J2HJ) has been most commonly adopted in previous calculations of this type. In all of our 
calculations, we estimated the contribution to the potential energy attributed to particles 
outside the main simulation cell by assuming that the pair correlation function g(r) equals 
one outside the cell. 
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Table U summarizes our results for the energetics of bulk 4 He, at different temperatures. 
Shown in parentheses are the corresponding results from Ref. The estimates are in quan- 
titative agreement, taking into account the statistical uncertainties of the two calculations. 
Amusingly, our statistical errors on the kinetic energy are not much smaller than those of 
Ref. 0, in spite of the fact that our calculation benefits of two more decades of advances in 
computing hardware. This is not too surprising, however, as the calculation of the kinetic 
energy, especially based on the estimator ()23|) . is known 13 to be the place where the limita- 
tions are most evident of using less than optimal a high-temperature density matrix p Q , such 
as the one used in this work. Still, the results of our calculation seem altogether satisfac- 
tory, giving us confidence that our PIMC implementation, including the new permutation 
sampling engine, performs correctly. 

B. Single-particle diffusion in imaginary time 

We observe excellent agreement between our results and those of Ref. Q for structural 
properties, such as the pair correlation function (an example is given in Fig. 1); however, 
effects of quantum statistics on these quantities are small,— and therefore their computation 
does not provide a particularly significant test of an algorithm to simulate indistinguishable 
quantum particles. 

More telling are measures of the diffusion of particles in imaginary time. Fig. 2 shows 
results for the quantity -D(t), defined as 



where (...) stands for statistical average. The two curves shown in the figures represent values 
of D computed by PIMC for bulk 4 He (solid line), as well as for a system of distinguishable 
4 He atoms, both at a temperature T=2 K. While in the first case 4 He atoms are treated 
as bosons, and therefore permutations are included, in the second case no permutations of 
particles are allowed. 

Obviously, because in the latter case one must have r(/5) = r(0), i.e., single-particle paths 
must close onto themselves, it must be D(j3) = 0. On the other hand, if permutations are 
allowed, then single-particle paths can become entangled, and D(f3) may take on a finite 
value. Moreover, the value of D is greater, at all imaginary times, in the case of Bose 





(25) 



6Ar 
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statistics; this is fairly intuitive, as the fact that particles are indistinguishable enhances the 
degree of delocalization of each individual particle. 

C. Superfluid density 

We have also computed the 4 He superfluid density ps, using the well-known "winding 
number" estimator.— At the lowest temperature considered in this work, namely T— 1.1 765 
K, our result is 1.02 ± 0.10, which is in agreement with experiment and with the PIMC 
result of Ref. Q. We obtained this result with a number of slices L = 544. It should also 
be mentioned that it appears possible to obtain a reasonably accurate estimate of ps using 
considerably fewer imaginary time slices (the result obtained with L = 136 is indistinguish- 
able, within statistical uncertainties, from the one quoted above), and that reducing L also 
causes a significant reduction of the statistical error on p$- In general, however, if L becomes 
relatively large, namely of the order of a few hundred, lengthy simulations are required in 
order to reduce statistical error to an acceptable size (e.g., 0.05 or less). This problem seems 
common to other PIMC implementations as well, and it is not clear to us to what extent it 
may signal an inefficiency of our sampling method. 

D. Statistics of Permutations and Permutation Cycles 

In order to characterize the performance of the permutation sampling algorithm, one 
may also look at quantities easily accessible in a simulation, which may not directly relate to 
anything measurable but provide a possible baseline for comparison of different algorithms. 
Table 2 provides statistics of permutation acceptance for a PIMC simulation of 64 4 He 
atoms at T— 1.1 765 K (the lowest temperature considered here). The total number of 
permutations attempted in this run is 4.5 x 10 5 , and the fraction of accepted permutations 
(of any cycle length) is approximately 0.4%. Permutations were sampled over an imaginary 
time interval of length 0.1 K _1 . 

As one can see from the second column, 2-particle permutations are sampled overwhelm- 
ingly more than others; however, the rate of acceptance of attempted permutations is essen- 
tially constant, independent of n. This is found to be the case at all temperatures considered 
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in this study. One may think that it would be advantageous to increase the rate at which 
permutations of more than two particles are sampled, since they presumably enhance the 
diffusion of the random walk through path space. Indeed, it is straightforward to gener- 
alize our sampling algorithm, so that permutations including more than two particles will 
be sampled more often. In practice, however, we found that including pair permutations is 
beneficial, in that it leads to a greater overall rate of acceptance of attempted permutations. 
How much of this is problem- or algorithm-dependent is difficult to say. Both quantities 
shown in Table 2 are rapidly decreasing functions of the temperature, as expected. 
Although they are not sampled directly, permutation cycles involving large numbers of par- 
ticles can and do occur, as a result of sampling many permutations involving few particles. 
Fig. 3 shows a histogram of probability for a particle to be part of a permutation cycle of 
length n (i.e., involving n particles) in a PIMC simulation carried out with the methodology 
illustrated above, at three different low temperatures below the A-transition. Although the 
data are somewhat noisy, these results are in quantitative agreement with those of Ceperley^ 
for the same system, using the CP methodology. As the temperature is lowered, the proba- 
bility that a particle will belong to a cycle of length n becomes independent of n. 



IV. CONCLUSIONS 



A new algorithm to perform the sampling of permutations of indistinguishable particle in 
Path Integral Monte Carlo simulations was introduced. This procedure is similar, in spirit, 
to existing methods, but differs in some important aspects, and may have some advantages. 
We have tested it by performing a PIMC simulation of liquid 4 He at low temperature, 
in the superfuid regime. Aside from the permutation sampling scheme, the rest of the 
PIMC methodology utilized here is not optimized for 4 He calculations. In particular, it is 
worth repeating that much better options exist for the high-temperature density matrix, 
which can drastically reduce the number of time slices needed for convergence. Still, the 
calculation proves quite feasible with currently available, moderately powerful workstations. 
It should also be noted that, while the use of a more accurate high-temperature density 
matrix (specifically, the pair-product approximation (j^J) greatly enhances the efficiency of 
calculations for a highly quantal, hard-sphere- like system such as helium, for other condensed 
systems such as molecular hydrogen, which feature a lesser degree of zero-point motion, or 
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Coulomb system, for which the interaction potential is considerably less "stiff", the high- 
temperature density matrix utilized here and Eq. (J5J) may be of comparable efficiency 
For comparison with existing calculations we limited the size of the system studied to N=6A 
particles, but it should be mentioned that simulations of systems with as many as four times 
more particles are also possible, with a reasonable amount of computer time (of the order 
of a month per thermodynamic point). 

We have attempted to furnish as much quantitative information as possible, that may help 
assess the relative efficiency of the permutation scheme proposed here against existing ones. 
Obviously, a direct comparison of results provided by implementations only differing by the 
permutation scheme adopted, is also desirable. It is our hope that such a comparison will 
be soon carried out. 
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TABLES 



TABLE I: Kinetic and potential energy per atom of condensed bulk 4 He, computed by PIMC at 
different thermodynamic conditions. Third column shows the number of imaginary time slices 
utilized in the calculation. The simulated system comprises 64 particles. Results in parentheses 
are from Ref. 2. 



T(K) 


p (A" 3 ) 


L 


kinetic (K) 


potential (K) 


1.1765 0.02182 


544 


14.123 ±0.028 


-21.3127 ±0.0025 








(14.17 ±0.08) 


(-21.35 ±0.04) 


1.379 


0.02182 466 


14.201 ±0.024 


-21.3195 ±0.0029 








(14.23 ± 0.08) 


(-21.35 ±0.08) 


1.600 


0.02183 400 


14.334 ± 0.034 


-21.3435 ± 0.0037 








(14.40 ± 0.08) 


(-21.39 ±0.04) 


1.818 


0.02186 


352 


14.468 ± 0.059 


-21.3913 ±0.0060 








(14.71 ± 0.08) 


(-21.44 ±0.04) 


2.000 


0.02191 


320 


14.862 ± 0.071 


-21.4810 ± 0.0078 








(15.05 ±0.08) 


(-21.57 ±0.04) 


2.353 


0.02191 


272 


15.821 ±0.062 


-21.5791 ± 0.0072 








(15.75 ± 0.08) 


(-21.60 ±0.04) 
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TABLE II: Permutation statistics for a PIMC simulation of bulk liquid 4 He at T=1.1765 K and a 
density p = 0.02182 A~ 3 . The second column indicates the percentage of attempted permutations 
involving n particles, sampled as explained in the text. The third column yields the percentage of 
attempted permutations that are accepted. The total number of attempted permutations for this 
run is 4.5 x 10 5 . 



n 


% attempted 


% accepted 


1 


0.0 


0.0 


2 


82.2 


0.38 


3 


12.5 


0.46 


4 


3.5 


0.27 


5 


0.1 


0.54 


> 5 


1.7 


0.00 
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FIGURE CAPTIONS 



FIG. 1: Pair correlation function g(r) computed by PIMC for liquid 4 He at a density p=0. 02182 
A~ 3 and at a temperature T=1.1765 K. The simulated system comprises 216 particles. Statistical 
errors are smaller than the size of the symbols. The effects of quantum statistics are barely visible 
on the scale of the figure. 

FIG. 2: Imaginary time diffusion coefficient D(t), defined in Eq. (|25|). computed by PIMC for bulk 
4 He (solid line) and for an ensemble of distinguishable (i.e., no permutations allowed) 4 He atoms 
at the same density (dashed line). The temperature is T=2 K. Statistical errors on both curves 
are of the order of 10~ 3 or less. The number of particles in the system is iV=64. 

FIG. 3: Probability for a single particle to belong to a permutation cycle including n particles, 
in a PIMC simulation of 64 4 He atoms at saturated vapor pressure and at three different low 
temperatures. 
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